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Abstract 

We present a mathematical formulation of kinetic boundary con- 
ditions for Lattice Boltzmann schemes in terms of reflection, slip, and 
accommodation coefficients. It is analytically and numerically shown 
that, in the presence of a non-zero slip coefficient, the Lattice Boltz- 
mann flow develops a physical slip flow component at the wall. More- 
over, it is shown that the slip coefficient can be tuned in such a way 
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to recover quantitative agreement with analytical and experimental 
results up to second order in the Knudsen number. 

1 Introduction 

Over the last decade, discrete kinetic methods, and most notably the Lattice 
Boltzmann (LB), have known a significant growth for the simulation of a 
variety of complex flows pQ . One of the most valuable properties of the LB 
method is its flexibility JT] , and easy set up of boundary conditions for com- 
plex geometries. Such a flexibility stems from the fact that the LB dynamics 
proceeds along rectilinear trajectories, so that LB shares the computational 
and conceptual simplicity of particle methods. On the other hand, since the 
LB dynamics typically involves more dependent variables (discrete distribu- 
tions) than hydrodynamic fields, the mathematical formulation of boundary 
conditions is left with some ambiguity, and a systematic treatment of the 
subject is still lacking. As a result, the issue of boundary conditions has 
become one of the most active areas of recent LB research. This is especially 
true for the highly debated question of the applicability of LB methods to 
microfluidic applications |2j. As recently shown by Ansumali and Karlin, [3], 
much can be gained in patterning LB boundary conditions after the time- 
honored procedures developed in continuum kinetic theory ^j. Based on the 
same philosophy, in this work we present a general class of homogeneous and 
isotropic boundary conditions for lattice Boltzmann models living on regular 
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lattices. This class represents a lattice realization of most popular boundary 
conditions in continuum kinetic theory, that is diffuse boundary conditions 
with and without accommodation. In this work, special emphasis is placed 
on the issue of slip flow at the boundary, which plays a central role in mi- 
crofluidic applications Analytical expressions for the slip flow are derived 
for a broad class of boundary conditions and compared with numerical sim- 
ulations. Excellent agreement between numerical and analytical results is 
found over a wide range of parameters. 

2 Formulation of the boundary condition prob- 
lem 

For the sake of concreteness, we shall refer to the two-dimensional nine-speed 
D2Q9 model [TQ| , although the proposed analysis can be extended in full 
generality to any other discrete-speed model living on a regular lattice (for 
example, the three-dimensional D3Q19 scheme shall be used for comparison 
with numerical details). 

We begin by considering the lattice Boltzmann equation in the following 
form: 

f i (x + c i 6 t ,t + 5 t ) - fi&t) = --(£(£,*) ~ ft Q \p^)) + %F 9i (I) 
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where q (i — 0, 1, 8) is a discrete set of velocities: 

Co = (0, 0)c 

ci,c 2 ,c 3 , c 4 = (l,0)c, (0, l)c, (-l,0)c, (0, -l)c 
c 5 ,c 6 ,c 7 , c 8 = (1, l)c, (-1, l)c, (-1, -l)c, (1, -l)c. 

The external source must inject zero mass and pF units of momentum per 
unit volume and time. This results in the following constraints on the forcing 
coefficients gf. 

Y^9iCi = YldiCi- Ci = 1. 

i i 

These constraints can be satisfied by choosing g 1 = —g 3 and g 5 — g 8 — —g-j = 
—g e , thus leaving us with one free parameter, say g 5 . 

We wish to emphasize that since the forcing term is associated with external 
forces, the ^ • || operator in the Boltzmann equation, there is no reason why 
g 5 should take the same value in the bulk and at the boundary. As a result, 
g 5 can be treated as an additional degree of freedom to describe fluid-wall 
interactions. 

The amplitude of the hydrodynamic forcing is chosen in such a way that, if 
the boundary velocity is zero, it reproduces a Poiseuille flow with centerline 
speed u : 
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where H = N y 5 y is the channel height. Discrete space and time increments 
are S x = 5 y ,S t , with c = ^, and the equilibrium distribution is given by: 

Aeq), r, . (3 ' u) . 1 (Q " «) 2 1 U 2 

with w = 4/9, iui,2,3,4 = 1/9, 1^5,6,7,8 — 1/36 and = J2i w i c "i — c 2 /3 the 
lattice sound speed. We consider a system whose y-coordinates lie between 
[0, N y 6y], N y being the number of grid points with solid walls located at —\5 y 
(south wall), and N y 8 y + \8 y (north wall). In the following, we will refer to 
the north wall in order to study the boundary condition problem. A general 
way of formulating the boundary conditions is: 

fM^ZjMMM*) ( 4 ) 

i 

where the matrix JCjj is the discrete analogue of the boundary scattering 
kernel expressing the fluid- wall interactions. In the above y = x + &i runs over 
the surface of the wall and the indices i,j stand for incoming and outgoing 
velocities respectively (j = 4, 7, 8 and i = 2, 5, 6 for the specific case of 
D2Q9)). To guarantee conservation of mass and normal momentum, the 
following sum-rule applies: 

E^ = ! (5) 

3 



5 



and upon the assumption of stationarity of the fluid-wall interaction, we can 
drop the dependence on t and write: 

fj = S^,i/i- ( 6 ) 

i 

We now focus our attention on the isotropic homogeneous case. The most 
general form for JCj^ is: 

lC jii = lC jti (\ci-c j \,\ci-n\) (7) 



where n is the outward unit normal to the surface boundary. The dependence 
on the second argument is necessary to develop a non-zero slip component 
in the stream-wise direction. 

For the case of D2Q9 (same goes for D3Q19), the quantity |q • n\ can take 
only two values, related to the two only possible angles that can be formed 
between a generic incoming velocity and the normal n, in our case a,\ = |7r 
and «2 = 7T- Explicitly, 



^ fl(x, Ny5y + 8y) ^ 

f4(X,N y 5y + 5y) 
\ f 8 {X,N y 5y + 5y) 



' h{x - 5 x ,N y 5y) 
f 2 (x,N y 5 y ) 

f 6 (x + 6 X , Ny5y) j 



(8) 
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where 



fClL 



OL2 



7T,Q2 



JCn. 

2 



Ql 



/C|, Q1 /Cs^c*! A^ai 



(9) 



Under the assumption of Stokes-flow, it can be shown that in the limit of 
zero spacings, the above scattering kernel provides an analytical expression 
for the slip velocity that, in the continuum limit of small time and space 
increments, reads: 



Usiip = A K Kn 



dn 



+ B K Kn 



wall 



On 2 



(10) 



wall 



where 



B 



K. 



12<7i 



1 + /C Ql - /Ci„, / c. 



+ 



"1 ^2' ai ' 
1 — fcn,a! + £%,ai 



12^5 



(1 + /C - ,ai) \1+JC 

7T,C 

In the above, rj = is the relaxation time of the continuum Boltzmann 
equation (proportional to the Knudsen number) and n is the inward normal 
to the wall. The case of finite spacings can be treated exactly (see Appendix), 
but for sake of simplicity, here we report only the continuum case. 
We wish to emphasize that what we call here 'slip velocity' is the velocity at 
the nodes nearest to the wall ((x, —\S y ) and (x, N y 5 y + \5 y ))\ to compute the 
velocity at the wall, an extrapolation of the parabolic profile at ( 
required. 
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Our analytical formula provides the slip velocity dependence on the stresses 
at the wall, on the Knudsen number and on the matrix elements of the 
scattering kernel, as well as on the forcing weights along the tangential and 
diagonal directions, g\ and g 5 respectively We now proceed to a more direct 
physical interpretation in terms of reflection and accommodation coefficients. 

2.1 Slip-Reflection model (SR) 

The first example involves two parameters r, s, representing the probability 
for a particle to be bounced back and slipped forward, respectively. The 
boundary kernel takes the form (see [S]): 



K 



r s 
r + s 
s r 



(11) 



Obviously, the two parameters are not independent and must be chosen such 
that r + s = 1. In this case the slip velocity (fTUj) reads : 



u 



slip 



AKn 



du x 
dh 



+ BKn 



wall 



d 2 u 2 
On 2 



(12) 



wall 



where 



.4 



c 1 — r 



c s r 



B = -(l-4g 5 ). 



(13) 
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From the first expression, it is clear that the coefficient A is equivalent to 
Maxwell's first order slip velocity 



a 



-Kn 



a 



du x 
dn 



wall 



with a = 2r [7] the well known Maxwell accommodation coefficient. It 
should be noted that in the limit of pure bounce-back (r = 1) the leading 
term disappears, and one is left with a purely quadratic dependence on the 
Knudsen number. This quadratic term stems from the tangential populations 
fi and ^3; the independence of the coefficient B of r is due to the fact that 
the tangential populations are evolved according to the same LB dynamics 
as in the bulk. We also note that in the limit of pure slip r — > 0, the slip flow 
tends to diverge, as it must be since this limit corresponds to zero friction at 
the wall. This contrasts with previous results [Hj, which reported a finite slip 
length even at vanishingly small values of r. The explanation is that those 
results were not converged in time. 

It is instructive to compare (fT^j) with the corresponding analytical solution 
for the fully continuum case (also with a continuum velocity phase space), 
namely ^ 



Us 



l.U6Kn 



dn 



+ 0.907Kn 2 



wall 



d 2 u x 



dn 2 



wall 



First, we notice that there exists a value of r such that the leading term can 
be reproduced exactly, that is: r* ~ 0.603. The quadratic term can also be 



reproduced exactly by choosing g 5 = 0.175 and g\ = 0.15. Note that this 
choice does not affect the bulk behaviour where the Stokes equation is still 
valid. 

Experimental data on slip flow are sometimes interpreted by assuming a 
Stokes flow in the bulk, coupled to a second order boundary condition of the 
form (|12j). In fact, in this way, upon the integration of the Stokes equation 
with the second order boundary slip it can be shown that the mass flow rate 
Q c of the channel is given by: 

Qc = SQ P (14) 

with Q p the Poiseuille mass flow rate and S a dimensionless number called 
slip coefficient depending on the Knudsen number : 

5=1 + 6AKn + UBKn 2 . 

In recent experimental works [S] the slip coefficient has been calculated for 
Helium and Nitrogen and a good experimental fit has been obtained with a 
second order slip at the boundaries that gives A = 1.2 ±0.05, B = 0.23 ±0.1 
for Helium, and A = 1.3 ±0.05, B = 0.26 ±0.1 for Nitrogen. It is interesting 
to note that both sets of coefficients can be exactly reproduced by our model 
by choosing r ~ 0.59, #5 ~ 0.22, (71 ~ 0.06 and r ~ 0.59, #5 ~ 0.23, g\ ~ 
0.04 for Helium and Nitrogen respectively. More generally, being subject to 
the constraint #5 < 1/4, our model permits to reproduce second order slip 
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coefficients in the range < B < 3. 



2.2 Slip-Reflection- Accommodation model (SRA) 

The second example is a straightforward generalization of the previous model, 
which is characterized by a third parameter, a, related to the presence of 
wall-relaxation (accommodation) phenomena at the fluid-solid interface. By 
accommodation we imply that the energy of the incoming and outgoing par- 
ticles is not the same because the outgoing ones are re-injected into the bulk 
with equilibrium weights. The SRA boundary kernel reads: 



/ 



r + aW 2 aW 2 s + aW 2 

aWi r + s + aWi aW\ 
s + aW 2 aW 2 r + aW 2 



\ 



(15) 



with the normalization constraints: 

r + s + a = 1 

and W 2 = 1/6, W\ = 2/3. The presence of the weights W\ and W 2 is due to 
the fact that they are the discrete analogue of the perfect accommodation 
kernel, i.e. a uniform Maxwell distribution at wall temperature. The SRA 
slip velocity has exactly the same form as (|T3*j) with the plain replacement 



r + a/2. 
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As expected, this implies that the accommodation coefficient a results in a 
smaller slip flow and this shows that the SRA model is basically equivalent 
to the SR one. 

It is known from continuum kinetic theory that a single accommodation 
coefficient is not sufficient for the quantitative interpretation of experimen- 
tal data. To cope with this problem, nearly four decades ago Cercignani and 
Lampis proposed a generalization of Maxwell's diffuse-boundary model which 
includes two accommodation coefficients, normal and tangential to the walls 
jH]. The lattice analogue of the Cercignani-Lampis kernel is readily com- 
puted. However, the analysis shows that only the accommodation along the 
tangential direction plays a role, while the one along the normal direction 
disappears. This is a pathology of all lattice models where incoming and 
outgoing normal velocities have the same magnitude. Therefore, we shall 
not consider this model any further in this work but the use of the lattice 
Cercignani-Lampis kernel might bear an interest for the case of thermal LB 
schemes. 

3 Numerical validation 

We now present a numerical study aimed at validating our analytical results 
for the SR and SRA kernels. To this purpose, we have performed numerical 
simulations of a channel flow with the 3d lattice Boltzmann model with 19 
discrete speeds (D3Q19). First, we consider the case of equi-balanced exter- 
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nal forcing among all biased populations (the equivalent of g$ = g\ in the 
case of D2Q9). The channel has dimensions N x = 64, N y = 32, N z = 32, 
the Knudsen number is Kn = 0.08 and the forcing has been fixed so as to 
reproduce a center-channel velocity Uq = 0.03 in the limit of a Poiseuille flow 

In Figure E]we show the center-channel profile in the stream- wise direction 
for the RS case. By varying the free parameter r between 0.1 and 1 the pro- 
file "shifts" with no changes in its concavity, which is fixed by the external 
forcing term F and the viscosity u, both kept constant in the simulations. 
In Figure 121 we present a comparison between the slip velocity, as extracted 
from numerical data, and our analytical results for the discrete case (see Eq. 
(|29|) in the Appendix). An excellent agreement between numerical and an- 
alytical results in the range 0.1 < r < 1 is clearly observed. For sake of 
completeness we have also compared the mass flow rate normalized to the 
pure bounce back case with our analytical solution (inset of Figure |2J) and, 
as expected, also here an excellent agreement is found. 
In Figure El we show the effect of the accommodation parameter a in reduc- 
ing the slip flow. The parameter is chosen as a = 0.3 and the values of r, s 
are chosen accordingly. The slip velocity is the same as the case SR, only 
with a renormalized reflection coefficient r = r + a/2 as can be clearly seen 
in the inset. 

We have also performed a set of numerical simulations with different repar- 
titions of the external forcing term, in order to study the slip properties 
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as a function of the Knudsen number. In the case of SR kernels, fixing 
the bounce-back parameter to r = 0.59 and the center-channel velocity to 
Uq = 0.01 in the limit of Poiseuille flow, we have studied the slip velocity as 
a function of the Knudsen number Kn in the range < Kn < 0.8. We have 
chosen an unbalanced bipartition with a forcing coefficient equal to 0.23 for 
the incoming populations in order to reproduce (see Figure HJ) the experimen- 
tal results showed in jH] for the case of Helium; the numerical results confirm 
our analysis and are indistinguishable from experimental data up to second 
order terms in the Knudsen number. 

4 Conclusions 

Summarizing, we have presented a unified formulation for a broad class of ho- 
mogeneous and isotropic boundary conditions for lattice Boltzmann models 
living on regular lattices. Analytical expressions for the slip flow at the fluid- 
solid boundary have been derived and successfully compared with numerical 
simulations of three-dimensional channel flow. The main conclusion is as 
follows: by allowing a non-zero slip coefficient, the Lattice Boltzmann flow 
develops a slip flow component which can be matched exactly to analytical 
and experimental data up to second order in the Knudsen term, well inside 
the transition regime (0 < Kn < 0.8). This means that the lattice Boltz- 
mann scheme with kinetic boundary conditions can be used to predict slip 
flow at finite Knudsen numbers well beyond the strict hydrodynamic limit. 
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Of course, the extension of the present model to more realistic situations, 
involving complex geometries and/or inhomogeneities ^21 G2] and thermal 
effects [Hj remains an open issue. 
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5 Appendix 

We impose stationary condition on the node (x, N y 5 y ) and, under the as- 
sumption of homogeneity, we drop the x dependence and write : 

*Wb(WA) + ^,aj2(N y S y ) + }C ha J 6 (N y S y ) - ^Fg 5 = f 7 (N y 5 y ) 
< K hai h{N y 8 y ) + IC^ a J 2 (N y 5 y ) + K^U{N y 8 y ) + $Fg 5 = f 8 (N y S y ) 
(/i " h)(N y 5 y ) - Uf q \p,u) - ft q \ P ,u)) = 2Ft$ 9i . 

(16) 

We can now write for the velocity in the x-direction at the height (N y S y ): 

\pUsii P = (2g 1 r + 2g 5 )F— + c(l - K^ ai + fcz tai )C(N y 6 y ) ( 17 ) 

with 

Mx,y)-f 6 (x,y) = {(y). (18) 
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Using the equations for f 5 and / 6 in the bulk, and under the stationarity 
assumption, we obtain: 

f 5 (x + c5 t , y + cS t ) - fc(x, y) = -±(f 5 (x, y) - ft v \p, u)) + 5 ^Fg 5 

f 6 (x - cS t ,y + cS t ) - fo(x,y) = -±(f 6 (x,y) - f { & eq \p,u)) - ^Fg 5 . 

(19) 

Dropping the x dependence (homogeneity), upon the definition (jTgj). we have 



C(y) - C(y - c5 t ) = -k(y - <*) + pUx[ l c6t) + 2F g ,% (20) 

ore cr 



In the limit S x , S y , St — » we obtain an O.D.E. which can be solved exactly 
However, by considering finite spacings such that S x — S y — c — 1, we have: 



C(J) - CO' - 1) = --CU - 1) + PUx{ l l) + 2Fg 5 l<j<N y (21) 
t 6r 



where j stands runs over the channel height. This finite difference equation 
can be solved exactly. The solution can be divided in two terms: the homo- 
geneous term C an d a particular term ( part (j). For the homogeneous 
term we have: 

c hom u) = (i - iy (22) 

while for the particular term we use the method of variation of constants: 
= (l~l) J ' l3 j:(l-ly\^+2Fg 5 }. (23) 
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Summing the two contributions (|23|) and (|22|) we have an explicit form for 
the general solution of (|2"T]) 

where C is a constant to be fixed based upon boundary conditions. Since 
C(N y ) enters in (|T7jl. we have: 

and under the assumption of low-Knudsen numbers (N y /r » 1), C{N y ) is 
well approximated by ( part (N y ), since 

(1 _ i)JV s o. As a result: 

CW^l-i)""" 1 !'^^)" 1 ^^^, (26) 

which substituted in (fTTj) . and setting p = 1, returns the following expression 
for the velocity at the height (N y ): 

1 / wXy- 1 ^- 1 f l\- k u (k) 

u slip = (6T gi +6g 5 )F+~(l-^ ai +}C hai )[(l - -J E (i " -) rft+WFgs)}. 

(27) 
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Next we impose a low-Reynolds regime, by specifying the velocity field as a 
(symmetric) parabolic profile: 



Ux(j) = a J 2 + b 3 + u *np 0<j<N y 
b = —aN y since u x (N y ) = u x (0) 

since vdjjU x (j) = —F 



(28) 



a = -f 



and we obtain for the slip velocity: 



u 



slip 



where: 



(12r 9l + 12g 5 )F /l - /C* >ai + /C* 



{^7T,ai - A. | ,ai + -I- J \A-7r,ai - A,| ,ai + ±/ 



1 - 

rV r 



^ N Ny—l Ny-l 



E *n 1 - 



fc=0 



/9 = - 1 - - 



]_ N Ny—l Ny-1 



fc=0 v 



(30) 
(31) 
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FIGURE 1 

Velocity profiles with SR kernels. We plot the center-channel profile in 
the stream-wise direction as a function of the normalized distance from the 
wall Unor-m — y/H, being H the channel height. Different values of the pa- 
rameters (r, s) are considered. From top to bottom we plot the following 
cases: (0.4, 0.6), (0.5, 0.5), (0.7, 0.3), (0.9, 0.1). In all simulations the Knud- 
sen number is Kn = 0.08 and the forcing term has been fixed to reproduce 
a center channel velocity uq = 0.03 in the limit of a Poiseuille flow (J2j). 

FIGURE 2 

Slip velocity for the SR kernels. We plot u s u p vs r for different values of 
the bounce back parameter r. The data of numerical simulations (boxes) 
are compared with the analytical estimate (|29|) (dotted line). The agreement 
is excellent all over the range 0.1 < r < 1. Inset: numerical data ((boxes) 
representing the mass flow rate normalized to its pure bounce-back value 
(Qnorm) are compared with our analytical estimate (dotted line). 

FIGURE 3 

Velocity profiles with the SRA kernels. We plot the center channel profile in 
the stream-wise direction as a function of the normalized distance from the 
wall y n0 rm = y/H, being H the channel height. Different values of the param- 
eters (a, r, s) are considered. We analyze the following plots: (0.0, 0.5, 0.5) 
(top plot) and (0.3, 0.5, 0.2) (bottom plot). Inset: analysis for the slip veloc- 
ity in the channel in the case (a, r, s) with a = 0.3. The data of numerical 
simulations (boxes) are compared with the analytical estimate (J2~9~j) (dotted 
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line) with a normalized reflection coefficient r = r + a/2. The agreement is 



excellent all over the range 0.2 < r < 0.7. 



FIGURE 4 



Analysis for the slip coefficient (fT4"j) in the channel as a function of the Knud- 
sen number. The data of numerical simulations (circles) with a properly 
unbalanced repartition of the external forcing are compared with the experi- 
mental results (triangles) given in jHj and with the analytical expression (|14|) 
with A = 1.2 and B = 0.23 (dotted line). The linear full accommodation case 



as a function of the Knudsen number is studied with the same repartition of 
the external forcing as before. We compare the numerical data (circles) with 
our analytical expression (dotted line) including, as before, terms up to the 
second order in the Knudsen number. In all the simulations we have used a 
SR kernel with a bounce back parameter r = 0.59. 




wall 



) is also plotted for a comparison. Inset: the slip velocity 
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Figure 2: 
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